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Abstract. We study the out-of-eqidlibrium evolution of a strongly interacting 
quantum spin chain which is mapped on a system of hard rods that are coherently 
deposited on and removed from a lattice. We show that this closed quantum system 
approaches an equilibrium steady state which strongly resembles a microcanonical 
ensemble of classical hard rods. Starting from the fully coherent evolution equation 
we derive a Master equation for the evolution of the number of hard rods on the 
lattice. This equation does not only capture properties of the equilibriirai state but 
also describes the d5mamical non-equilibrium evolution into it for the majority of 
initial conditions. We analyze this in detail for hard rods of varying size. 
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A central objective of non-equilibrium statistical mechanics is the description of 
the time evolution of thermodynamic properties of macroscopic systems that are 
prepared far from equilibrium. In principle one should think that it is always 
possible to derive an effective equation of motion, e.g. a Fokker-Planck equation, 
that governs the non-equilibrium behavior of thermodynamic quantities, from the 
underlying microscopies. In general this is, however, a formidable task and even 
more so in isolated systems where there is no apparent system-bath decomposition 
or an immediately manifest hierarchy in the coupling strengths between two or more 
subsystems. Nevertheless, these isolated systems often exhibit a generic relaxation 
dynamics, in the sense that for a subset of observables the non-equilibrium evolution 
leads to a steady state which is well-described by a thermodynamic ensemble. 

For classical systems the underlying microscopic Hamilton equations of motion 
are non-linear in the dynamical variables. Generically, this leads to chaos such 
that the system ultimately explores the entire phase space. The occurrence of 
thermalization is then made plausible by invoking the Ergodic Hypothesis [1]. In 
contrast, the microscopic equation describing the dynamics of a closed quantum 
system is a linear wave equation for a complex valued wavefunction and observables 
are represented by Hermitian operators. Due to these mathematical differences it 
is not obvious which equivalent arguments would lead to thermalization in this 
quantum context. 

Given this, a first approach to study thermalization of quantum systems would 
involve a direct calculation of their microscopic dynamical evolution. However, 
determining the exact quantum dynamics of an interacting many-body system is 
a formidable task. The microscopic approach of solving Schrodinger's equation 
is usually restricted to very small systems, since the exponential growth of the 
Hilbert space with the number of particles also entails a growth in the mathematical 
complexity of the problem that quickly becomes unmanageable. Circumventing this 
problem requires sophisticated numerical tools and clever mathematical approaches. 
It is fair to say that our understanding of the equilibration and thermalization 
of closed quantum systems is still far less developed than for their classical 
counterparts. 

Recently, there is an ever growing interest in addressing the question of when 
and how these systems relax towards an equilibrium state when prepared in a non- 
stationary initial state [2-16]. This development is, on the one hand, sparked by 
the availability of state of the art experiments - especially in the field of ultracold 
gases - where the investigated systems are extremely well isolated from the thermal 
environment and thus constitute almost ideal realizations of closed systems [17-21]. 
These experiments do not only constitute ideal test-beds for new theoretical ideas and 
concepts but they also point to the fact that there is growing interest for a quantum 
theory of the dynamics of closed systems to adequately interpret experimental 
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observations made in these thermally isolated setups. 

On the other hand, recent successes in the theory of equilibration and 
thermalization of closed quantum systems have fueled the interest in this 
topic. A very insightful idea was formulated with the Eigenstate Thermalization 
Hypothesis (ETH) that links thermalization to the spectral properties of the quantum 
Hamiltonian [22, 23]. In essence, ETH states that in a thermalizing system any 
many-body eigenstate contains a thermal state, i.e., that the expectation value of an 
observable O calculated in an eigenstate with energy e coincides with the thermal 
average taken in the microcanonical ensemble at same energy. This conjecture has 
numerically been verified for a number of model systems [24-26]. 

Despite its success the ETH does not provide an answer to how relaxation takes 
place and does not provide any information on the non-equilibrium evolution which 
eventually drives the system into a steady state. Here we take a closer look at the 
equilibration of a closed, interacting many-body quantum system in the time domain, 
expanding on a recent work [16]. To this end, we consider a simple yet generic one- 
dimensional spin model, which belongs to the class of Ising models in a transverse 
field and can be mapped on a system of quantum hard rods that are coherently 
removed from/deposited on a lattice. For this system we show that it is possible to 
derive analytically a Master equation which captures the non-equilibrium dynamics 
of the number of hard rods. We undertake a comparison of this effective description 
with numerically exact simulations of the fully quantum problem and find excellent 
agreement. In this way, the results of this work also contribute to recent efforts to 
derive effective diffusion equations that describe the non-equilibrium dynamics in 
closed quantum systems [16,27,28]. 

This paper is organized as follows: In section 2 we introduce a spin Hamiltonian 
that represents a quantum version of a lattice gas of hard rods. Subsequently, 
in section 3, we introduce a specific graphical representation of the Hilbert space 
of the quantum model - the configuration network - to conveniently illustrate the 
dynamics of this system. Using the insights from a statistical analysis of the network, 
we present the derivation of a Master equation that describes the dynamics of the 
number distribution function of hard rods in section 4. In particular, we show that 
the rate coefficients of this Master equations can be calculated analytically. Moreover, 
we numerically investigate the time evolution and steady state distribution of the 
density of hard rods. Following this analysis, we compare in section 5 the predictions 
of the Master equation with the results obtained by numerically solving the many- 
body Schrodinger equation for our model system and find excellent agreement. 
Conclusions and an outlook are provided in section 6. Unless stated otherwise, we 
set /i = 1 in this work. 
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2. The system 

The quantum system in which we are interested here is a chain of L equally spaced, 
interacting spin- 1/2 particles at nearest neighbor distance a in a transverse magnetic 
field of strength Q. Its Hamiltonian is given by 



where, crj = (|t) {l\ + I4-) {tDj is a Pauli spin matrix effectuating spin flips of the j-th 
spin and uj = (|t) (t|)j is a projector on the up-state (excited state) of the j-th spin. 
The spin-spin interaction of strength V depends on the spin-state as well as on the 
separation of the spins: two spins interact (i) if both of them are simultaneously in 
the spin-up state and (ii) if they are separated by a distance that does not exceed a 
critical distance A, which we will refer to as blockade radius. We use this terminology 
because in this work we focus specifically on the regime V/Q, ^ oo. This means that 
it is energetically forbidden to have two or more up-spins located within the critical 
radius Tc = A a, i.e. an up-spin blocks the excitation of spins to the up-state in its 
vicinity. The model defined by the Hamiltonian (1) is in fact of practical relevance, 
as it has been shown to capture the quantum dynamics of strongly interacting, laser- 
driven Rydberg atoms trapped in an optical lattice [29,30]. For more details we refer 
the reader to references [31-34]. 

The Hilbert space of the system is spanned by all classical spin configurations 
compatible with the blockade condition. In figure 1(a) we provide examples of 
permitted spin configurations for the cases A = 1 and A = 2. Instead of using spins 
we employ a description of the basis states in terms of hard rods. Each configuration 
of spins can be uniquely mapped into and arrangement of hard rods of length A-M, as 
illustrated in figure 1(a). A spin-flip from the down- to the up-state corresponds to the 
deposition a hard rod, while the opposite process removes a hard rod from the lattice. 
Our aim is to analyze the non-equilibrium dynamics of the hard rods effectuated 
by i^spin after a quench. In particular we seek an effective evolution equation for 
the probability density Pn{t), which describes the probability of having the lattice 
occupied by n hard rods at time t. The first step towards such equation is to gain 
an understanding of the Hilbert space structure. Here it is convenient to employ a 
graphical representation of the Hilbert space as a configuration network. We will see 
that the structure of the effective evolution equation for the pn will crucially depend 
on the statistical properties of the network. In the following, we will construct the 
network and analyze its properties in detail. 
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Figure 1. (a) Representation of configurations in terms of spins (arrows) and the hard 
rods (grey rectangles). Each up spin inhibits further excitations on A neighboring 
sites on both its left and right. This blockade effect can equivalently be illustrated 
by hard rods on the lattice, each occupying A + 1 sites, and spin flips correspond 
to placing/ removing hard rods on the lattice, (b) An example of the configuration 
network for L = 8 and A = 1. Each node represents a classical arrangement of n 
hard rods |?iC„) (see text). The links between nodes connect configurations that can 
be converted into one another by the deposition/ removal of a hard rod. Since we 
are interested in the dynamics in number space we do not resolve individual states 
but consider the time evolution of the operator P„ that projects on all basis states 
contained in the n-th column. The equation of motion for (P„) depends directly on 
the structure of the graph, specifically on the abundance of loops (green), reflections 
(orange) and transmissions (red). 

3. Configuration network 

3.1. Structure 

To construct the configuration network we introduce the basis states |nC„). Each 
of these states represents a specific classical arrangement C„ of n hard rods. These 
states satisfy the completeness X]nc„ I'^^n) ('^Cnl = ^ and orthonormality relation 
(nCnlmlCm) = 5nmSc„K.m- interpret each of the microstates |nC„) as a node of 
a network. By grouping configurations containing the same number of hard rods 
into columns we obtain the network structure depicted in figure 1(b). The time- 
evolution of the system can then be imagined as a temporal change in the occupation 
of these nodes. Dynamics is introduced through the Hamiltonian Hq which leads 
to transitions between microstates that are represented as edges of the network. 
Since Hq causes only single spin flips, nodes in neighboring columns are only linked 
directly, if their corresponding microstates can be converted into one another by the 
removal/ deposition of one hard rod. For example, setting A = 1 the state | Itiitiiti) 
in the n = 3 column is directly linked with ||tiiiiiti)/ but not with liitiiiti) iri 
the n = 2 column. In the following we will analyze in detail the properties of the 



Equilibration of quantum hard rods in one dimension 



6 



conjfiguration network. 



3.2. Properties of the configuration network 

The most basic properties that define the structure of our configuration network 
are the number of columns and the number of nodes within each column. Fixing 
the length of the system to L sites, and applying periodic boundary condition, the 
maximum number of hard rods that each occupy A + 1 sites (i.e., the blockade radius 
is A), which can be placed on the lattice is [L/(A + 1)J, where [x\ denotes the closest 
integer smaller or equal to x. The index counting the number of hard rods can thus 
take the values n = 0,1,...,LL/(A + 1)J. The number of microstates z/„ contained 
in the n-th column is given by the number of ways in which one can distribute n 
indistinguishable hard rods of length A + 1 over L lattice sites. This is a standard 
combinatorial problem with solution 

_ L{L-1-Xn)\ 
''"~n!(L-(A + l)n)!- ^ ^ 

Having determined the properties of the "backbone" of the network we now turn 
to assessing the linkage of the nodes. In particular, we calculate the mean number of 
different possibilities r„^n±i to go from a state with n hard rods to one in the adjacent 
columns. This quantity can be expressed as T„^„±i = c„ „-ti/z/n, where c„ „-|-i denotes 
the total number of links between columns n and n± 1 and hence c„±i,n = c„,„±i. (Note 
that this symmetry does not hold for the T„^m.) Moreover, we know that T„^„_i = n, 
as in a configuration of n hard rods there are n possibilities for removing one hard 
rod reaching a state with n — 1 hard rods. Using these relations the total number of 
links between two columns evaluates to 

L{L-1-Xn)\ 



(n-l)!(L-(A + l)n)!' 

L{L - 1 - \{n + 
~ n!(L-(A+l)(n + l))!' 

(3) 



which can then be used to calculate T„_^„+i. 

Let us continue by analyzing the second order processes shown in figure 1(b), i.e. 
loops, reflections and transmissions, as they will enter in the derivation of the Master 
equation. Selecting a node from the network, the number of loop transitions from 
that node equals the number of links that this node has with other nodes. Therefore, 
the mean number of loop transitions from a state with n hard rods is given by 



-^/o'o'p — + (4) 

Reflections connect two configurations that contain the same number of hard rods but 
differ in the position of exactly one hard rod [cf. figure 2 (a,b)]. Two configurations 
that are randomly selected from one column will typically differ by the positioning of 
more than one hard rod and are thus not connected by a reflection. If two microstates 
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Figure 2. (a,b) Graphical illustration of the reflection processes. Two configurations 
can at most be connected by two reflection pathways (a). However, since rods must 
not overlap, the path containing an additional hard rod during the intermediate step 
is often forbidden (b). (c) For the numerical study we fix the physical system length 
I and the critical radius r^- The parameter A is varied by increasing the number of 
lattice sites. 



happen to be connected there are at most two paths as illustrated in figure 2: a 
deexcitation followed by an excitation or vice versa. For n ^ 1 (at high density) 
there is often even only one path available [figure 2(b)]. Similar considerations can 
also be made for transmission diagrams, i.e., the average number of paths connecting 
two microstates containing n and n±2 hard rods is ~ 1 . In fact the mean number of 

reflections (transmissions) A'^^^^'^^ (^traL) between two randomly selected states can be 
calculated analytically: 

f.T{n) >n(^7i+l— 5>n l)'^n+l ^n— 1— 5>n(^n— 1— >n l)'^n— 1 

~ l^nil^n - 1) l^nil^n " 1) ' 

njin) !>n+l^n+l— !>n ~l" ^n— 2— s>n— l^n— 1— !>n 

-''^ trans 



(5) 



To illustrate that loop transitions are far more abundant than reflections and 
transmissions we present some numerical examples in table 1. Here we compare 

^/oop' ^rlfi ^traL ^ number of lattice and hard rod sizes. This leads to two 
observations. First, the relative weight of loop transitions largely increases with 
increasing systems size L, as due to the larger dimension of the Hilbert space each 
state can have more connections to other configurations. Second, the probability 
that two states within a column of the network are connected by a reflection is 
vanishingly small in the "bulk" of the network (1 ^ n ^ [L/(A + 1)J). The same 
is also true for transmissions between columns n and n ± 2. Note, that near the 
boundaries of the network, i.e. columns close to the maximum /minimum n-value. 
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280.1 


0.026 
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93.39 


1.64 X 10-27 


3.18 X 10-27 
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244.78 


0.025 


2 




10 


68 


3.11 X 10-^3 


1.73 X 10-^2 



Table 1. Average number of connection between nodes established by the three types 
of second order processes depicted in figure 1(b). We have chosen a relatively small 
lattice size (i = 12) and a large one (L = 300) for illustration. 



the condition N{^^^ » -^trans) ^refi l^ss ^^11 satisjfied. However, this concerns only 
an exponentially small subset of states forming the configuration network. These two 
observations on the statistics of the configuration network are of central importance 
in the derivation of the effective Master equation for Pn{t) that we are going to present 
in the following section. 

4. Master Equation 

In this section we derive an equation of motion which describes the evolution of 
number of hard rods on the lattice as function of time. This means that we are 
not interested in the actual population of individual nodes of the graph shown in 
figure 1(b) but rather in the probability p„(t) of the system to reside in a specific 
column n at time t. We will see that this eventually leads to a Master equation 
that has a steady state Pn{t oo) which is proportional to the number of classical 
configurations contained in a specific column. This strongly suggests that this steady 
state corresponds to a microcanonical equilibrium state in which all arrangements of 
hard rods occur with equal probability. 

4.1. From the von-Neumann equation to the Master equation 

The probability p„ is defined as Pn = Tr p{t)Pn, where p{t) is the density matrix of 
the system and Pn = J2cn 1^^") is a projector which projects onto the subspace 
spanned by all microstates contained in the n-th column of the configuration network 
[see figure 1(b)]. Throughout this work we are interested in a situation where the 
initial state of the system p(0) contains a fixed number of hard rods, i.e. [Pn, p(0)] — 0. 

To begin with the formal derivation of the Master equation, let us momentarily 
return to the case of finite V and transform the Hamiltonian (1) into the interaction 
picture with respect to Hy. As shown in Appendix A one can then write the evolution 
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equation for p„ = (Pn) — Tr p{t)Pn as 

dt{Pn)t= - [ dsTT{PnHi{t)Hj{s)p{s) + Hi{s)Hjit)Pnp{s) 
Jo 

- Hi{s)PnHi{t)p{s) - Hi{t)PnHi{s)p{s)}, (6) 

with 

Hi{t)^ J2 {nC^\Hn\m}C^)e-'^''->'rr.---^r.)t\nCn){mJCm\, (7) 

where ujrvCr, — {'nCn\Hv\nCn) denotes the conjfiguration energy of the spin 
configuration |nC„). In order to evaluate the integral in the still exact evolution 
(6) we make the replacement p{s) — )■ p{t), which effectively amounts to a second 
order approximation. At this point the quality of this approximation is not clear 
but we will provide a numerical justification of this step a posteriori. Returning to 
the ideal blockade regime {V /Vt — )■ oo) the fact that each allowed spin configuration 
has the same configuration energy cu^^ = 0, removes the exponential factor of the 
interaction picture Hamiltonian making it time-independent {Hi{t) — )■ Hq). Thus, 
the time-integration in (6) can be carried out and simply amounts to a multiplication 
by t. The steps result in the second order equation 

dt{Pn)t = -tTr [PnHl + HlPr, - 2HnPnHn] p{t), (8) 

that depends on products of the form 

{nCn\Ha\mlCm) {mlC^\H^\pCp) = (nC„| Hj^ \pCp) . (9) 

niKm 

These matrix elements are the second order processes discussed in the previous 
section. Diagonal elements correspond to loops and off-diagonal ones to either 
reflections or transmission. From the analysis of the graph we know that loop 
transitions are far dominant and that Hq is thus approximately diagonal, hence 

{nCn\ \pCp) ^ {nCn\ |nC„) SnpSc„Cp- 

The final approximation consists of neglecting the variation of the matrix elements 
of Hq^ with Cn within a given column of the configuration network, by replacing the 
matrix element with its average taken over all microstates within a column 

{nCn\ Hi \nCn) ^ O^TV/^jp = ^^^^^^^ ^ Tr^^n-^) . (10) 

Inserting these approximations into (8) and writing Pn{t) = T^^PnPit) = 
Tr \nCn) {nCn \ p{t) we arrive at the desired Master equation for the dynamics of 
the probability distribution p„ of finding n hard rods in the system at time t 

it)] 

- 2nH [Tn^n-l + Tn^n+l] Pn{t). (11) 

This Master equation has a stationary solution despite the explicit appearance of the 
time variable on the right hand side. This can be seen by transforming to a new time 
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Figure 3. Time evolution of the distribution function p„ (t) as a function of the hard 
rod density nX/L for A = 1 (nearest neighbor blockade). The system considered has 
a lattice length to blockade radius ratio l/vc — 120. From left to right p„ is shown 
at: nt = 0.02(red), 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8 and f7i = 2 (blue, corresponding to the 
steady state distribution). The initial distribution function used is po = dn,o- 



variable r = fit^ which removes the explicit time-dependence. The rate coefficients of 
the Master equation (11) obey detailed balance, i.e., p^Tn^n+i = Pn+i^n+i-^n/ where 
denotes the steady state distribution 

eq ^ Q2) 



y 



=0 



This distribution is proportional to the number i/„ of arrangements of n hard rods, 
suggesting that each node of the network is populated with equal probability, as one 
would expect from a "maximum entropy state". 

As shown in reference [16], (11) is a discretized Fokker-Planck equation in 
excitation number space. Its drift coefficient is proportional to T„^„ i — T„^„+i and 
its diffusion coefficient proportional to T,j^„_i + T„^„+i. The Master equation thus 
describes the excitation dynamics of the system as diffusion between the columns of 
the configuration network shown in figure 1(b). 



4.2. Time Evolution and Steady State of the Master Equation 

We will now study the time evolution of the distribution function pn by numerically 
solving (11). We perform the following analysis by fixing both the physical lattice 
length / and the critical radius Tc. The ratio //rc is then also fixed. We choose a value 
oil/vc = 120 in the numerical examples of this subsection, which means that at most 
120 hard rods can be placed on the lattice. What we change is the number of lattice 
sites occupied by a hard rod, i.e. A + 1, which is done by varying the total number L 
of lattice sites. This is illustrated in figure 2(c). We note that for A — )■ oo this procedure 
results in taking the continuum limit a/r^ — )■ 0. 
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Figure 4. Dependence of the time evolution and steady state properties on A. (a) Time- 
evolution of the average hard rod density {n)X/L starting from the configuration 
without hard rods for different numbers of blockaded sites: A = 1 (blue), A = 9 
(green) and A = 59 (red), (b) Steady state values of {n)X/L as function of A. (c) Steady 
state distribution p„ as function of the hard rod density nX/L for the three values of 
the number of blockaded sites used in (a), (d) Steady state value of the parameter Q 
as function of the number of blockaded sites. For more details, see text. 

Figure 3 shows the time evolution of p„ as a function of the hard rod density 
nX/L for A = 1 starting from the initial state with no hard rods, i.e., Pn(0) = 5„_o- The 
time evolution of the distribution function towards its steady state (indicated in blue) 
shows the characteristic features expected from a system following a Fokker-Planck 
evolution: p„ gets broader due to diffusion between the columns of the configuration 
network and drifts towards its equilibrium position, where the rate coefficients for 
further excitation (r„^„+i) and de-excitation (r„^„_i) become equal. 

Figure 4 summarizes the behavior of the solution of the Master equation for 
varying A. Panel (a) shows the time evolution of the mean hard rod density {n)\/L 
starting from an empty lattice. For short times the hard rod density exhibits a 
quadratic time dependence, {n)\/L oc t^. This behavior is induced by the explicit 
linear time dependence of the right hand side of the Master equation. For long 
times in)\/L saturates to its steady state value, which is reached faster the larger 
A. Furthermore, the hard rod density in the steady state gets larger as A increases 
[figure 4(b)]. Interestingly, however, the fluctuations around the mean density {n)\/L 
decrease with increasing A, which can be seen in figure 4(c), where the full steady 
state distribution function p„ is depicted for three different values of A. This behavior 
is further quantified in figure 4(d), where we show the steady state value of the ratio 
Q = {{n^) — {n)'^)/{n), which relates the width of the distribution function to its 
mean, as a function of A. The ratio asymptotically approaches zero with increasing 
A, i.e. increasing number of lattice sites occupied by a hard rod [cf. figure 4(d)]. This 
surprising feature is the result of an entropic effect which has recently been discussed 
in reference [33]. 




Figure 5. Comparison of the exact quantum d5mamics with the predictions of the 
Master equation (11) for different A. (a) Time evolution of the hard rod density and (b) 
the distribution functions at iU — 20. The vertical dashed lines indicate the maximum 
possible hard rod density for each A, which is given by A/(A + 1). From top to bottom, 
the blockade radius is increased as A = 1,2,3 while the ratio l/vc = 10 remains fixed. 
In all plots, the black curves are the solutions of the Master equation. The colored 
curves show the solution of the Schrodinger equation using three randomly chosen 
initial states | noCng ) with the same number uq of initial hard rods (from top to bottom 
no = 2, 2, 3). 



5. Comparison with the Exact Quantum Evolution 

After having discussed the main features of the time evolution and the steady state 
of the Master equation we will now compare its predictions to the exact quantum 
dynamics of the system. To this end we have numerically solved Schodinger's 
equation with the Hamiltonian (1) in the limit oiV/^l oo for up to L = 30 sites. The 
left column of figure 5 shows the numerically exact quantum evolution of the hard 
rod density together with the prediction of the Master equation (in black). The ratio 
of system length to critical radius is fixed to l/vc = 10 and we choose A = 1, 2, 3 from 
the top to the bottom panel. Note that the dimension of the Hilbert space increases 
from top to bottom. The differently colored curves in each panel show solutions to 
the Schrodinger equation starting from randomly chosen initial states | raoC„o ) with the 



Equilibration of quantum hard rods in one dimension 



13 



same number of hard rods no but different spin configurations C„„ . The initial number 
of hard rods is chosen such that the initial state lies in a region of the configuration 
network with large connectivity, i.e., where the statistical assumptions underlying 
the derivation of the Master equation are well met. The right column shows the 
corresponding probability distributions p„ at f2t = 20. For both {n)X/L and p„ the 
agreement between the results of the exact quantum calculation and the prediction 
of the Master equation is remarkably good. In particular, for long times the results 
of the full quantum calculation and solution of the Master equation only differ by 
roughly one per cent. For short times the quadratic time dependence of {n)X/L as 
well as its dependence on A are well reproduced (see insets). For longer times the full 
quantum solutions exhibit oscillations around the equilibrium value of {n)\/L. Being 
a simple rate equation our Master equation does not reproduce these. However, the 
quantum oscillations decrease with increasing dimension of the Hilbert space. For 
A = 1 this behavior was also reported in reference [16]. 

At long times we observe a small systematic offset of {n)X/L obtained from the 
exact numerics from the steady state values predicted by the Master equation. In 
the cases shown in figure 5(a), where the initial state contains fewer hard rods than 
the equilibrium state, the quantum calculations suggest a sightly lower value of the 
average hard rod density at long times. In contrast, for initial states with a higher 
hard rod density than the equilibrium value the quantum results lie sightly above 
the prediction of the Master equation. Due to this systematic dependence on the 
initial state, we attribute this small offset to the presence of memory effects in the 
quantum dynamics that were completely neglected in the derivation of the Master 
equation. Furthermore, the data of the full quantum calculation for A = 3 exhibit a 
very slow drift of {n)X/L towards the steady state of our rate equation. This effect 
is seen more clearly in figure 6, where we follow the time evolution of the hard rod 
density for A = 4 to much longer times. The observed shift might be indicative of 
a pre-equilibration process in the quantum system, in which {n)X/L quickly reaches 
a quasi-equilibrium state, which then very slowly equilibrates to the "true" steady 
state. However, the drift might also stem from a long wavelength oscillation present 
in the full quantum calculation due to the finite size of the system. Since the full 
quantum calculations are limited to small system sizes it is difficult to further explore 
this effect, which seems to be more pronounced with increasing A. 

Let us finally return to the observation made in figure 5(b) that distribution 
functions pn calculated fully quantum mechanically and using the Master equation 
agree very well at long times. In order to quantify the degree of agreement we use 
the following overlap measure [35], 



Here p'^^ denotes the steady state solution of the Master equation and pn is the 
equilibrium distribution obtained from the solution of the Schrodinger equation. 




L/iX+l) 



(13) 
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Figure 6. Long time quantum evolution of the hard rod density starting from a 
randomly chosen initial state with no = 3 hard rods for A = 4 and l/vc = 9. The 
evolution predicted by the Master equation is plotted in black. To highlight the shift, 
the time evolution from fit = 0.5 to fit — 20 is enlarged and shown in the left inset, 
and the evolution from fit = 80 to fit = 100 is enlarged and shown in the right inset. 



time averaged over a time interval At, i.e., p„ = drp„(r)/At. The distribution 
functions are identical when V = 1, and V = for distribution functions that are 
completely non-overlapping. For the three situations discussed in figure 5 we have 
selected 100 randomly chosen initial spin configurations for L = 10, A = 1 and 500 
randomly chosen initial spin configurations for L = 20, A = 2 and L = 30, A = 3, and 
have used them as initial states of the quantum evolution. The time average in order 
to compute p„ was taken over the interval At = [20/^2, 40/i7]. We have collected the 
results of these simulations in the histograms shown in figure 7. Here we see that 
for the vast majority of initial conditions V is close to one with only a few outliers. 
This indicates that indeed equilibration is largely independent of the initial state. In 
order to demonstrate that the used set of initial conditions is representative we have 
looked at the distribution of the number of hard rods contained in the initial states. 
As an example the inset in figure 7 shows this distribution for the parameters of the 
bottom panel. Comparing this to the microcanonical steady state distribution we 
see that indeed each number is represented with the correct weight. Finally, looking 
at the histograms shown in figure 7 as a function of the dimension of the Hilbert 
space we see that the system equilibrates better with increasing number of available 
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Figure 7. Histogram of the parameter V, defined in (13), for the three parameter sets 
used in figure 5. The numbers in the parentheses show the dimension of the Hilbert 
space and the number of initial states used to create the histogram, respectively. For 
the bottom panel we also show the distribution of the number of hard rods contained 
in the initial states as an inset. The black dashed line with blue cross marks in the 
inset shows the microcanonical steady state distribution of the number of hard rods 
for comparison. 



microstates, confirming similar observations made in [16] for the case A = 1. In the 
thermodynamic limit the Master equation (11) therefore indeed provides an excellent 
description of the non-equilibrium dynamics in particle number space. 

6. Conclusions and Outlook 

In this work we have investigated the non-equilibrium dynamics and equilibration 
of a strongly interacting closed quantum system of hard rods of varying size. We 
showed that the steady state and the non-equilibrium relaxation towards it are well 
captured by a Master equation whose numerical coefficients can be analytically 
calculated using the combinatorics of hard rods placed on a lattice. The dependence 
of the relaxation dynamics on the initial state has been systematically investigated 
through numerical simulations of the full quantum dynamics. It has been found 
that relaxation into the steady state is achieved for the majority of initial conditions. 



Equilibration of quantum hard rods in one dimension 



16 




Figure 8. Configuration network of a system with L = 6 sites, A = 1 and finite ratio 
V/ft. Different layers represent different interaction energies and adjacent layers are 
separated by V. The lowest layer corresponds to a network similar to the one shown 
in figure 1(b). Due to the structure of the Hamiltonian transitions take place only 
between adjacent layers and layers that are separated by 2V. The timescale for inter- 
layer-transitions is ^ fl^/V while intra-layer-relaxation takes place on a time ^ fi"^. 
The index j labels the number of adjacent nearest neighbor excitation pairs. 

Expectation values taken in the steady state are compatible with the assumption of 
an equal population of each microstate, suggesting that the system here is indeed 
well described by a microcanonical ensemble. 

To conclude, let us finally comment qualitatively on large but finite interaction 
V, i.e. where the system no longer can be described in terms of hard rods. In this 
more general case we can still employ a representation of the state space in terms 
of a graph, but as shown in figure 8 (for A = 1) a third dimension has to be added 
which accounts for the interaction energy V of different combinations. The lowest 
layer corresponds to the set of nodes in the V/Vl oo limit centered around zero 
energy with no contiguous excitations, cf. figure 1(b). All higher layers contain j 
pairs of adjacent excitations and hence are centered around the energy ej = j x V. 
Each of these layers can be thought of as an energy shell for which an intra-layer 
diffusive Fokker-Planck dynamics can be derived. In reference [4] it was numerically 
shown that the steady state population of these energy shells follows a Boltzmann 
law pj = p^. (X e''^^^, with an inverse temperature j3, provided that V :§> Vl. This 
can be understood as follows: Choosing an initial state in the lowest energy shell 
the diffusive dynamics establishes a 'maximum entropy state' within this shell on 
a timescale fi^. Transitions to higher energy shells take place on a slow timescale 
~ ^I'^/V. Single spin flips as effectuated by the Hamiltonian can, however, only 
couple energy shells with a distance of at most 2V, i.e. a single spin flip can at 
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most create two pairs of adjacent excitations. That means that the population of 
j-th energy shell can be estimated by {Q'^/V^y which approximates a Boltzmann 
distribution with temperature such that ^ —2{^}/V)\og{^l/V). It remains an 
intriguing question whether also in this more complex system the derivation of a 
Master equation is possible from first principles. 
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Appendix A. Integrated von-Neumann equation 

In order to keep the notation compact we abbreviate our basis states as jo;) = \nCn)- 
The Hamiltonian given in (1) can then be expressed in this basis as 

H = J2 {a\H\a) \a){a\ + {o^W) \a){P\ ^ Hy + Hq. (A.l) 

We proceed by transforming the Hamiltonian into an interaction picture with respect 
to Hy by applying a unitary transformation U = exp[— itify] where we have set 
h = 1. The interaction picture Hamiltonian which is given by Hi{t) — U^HU, takes 
the following form: 

H:{t) = ha,0e-'^^'-^-^' \a) {P\ , (A.2) 

where cua = («|i/y|«) and ha^/s = {a\HQ\l3) are the diagonal and off-diagonal entries 
of the Hamiltonian, respectively. The density matrix p{t) of the system evolves (in 
the interaction picture) according to the von-Neumann equation 

dtp{t) = -z[Hr{t),p{t)]. (A.3) 

We proceed by formally integrating the von-Neumann equation, finding that, 

p{t)^p{0)-z fds[Hi{s),p{s)]. (A.4) 
Jo 

Substituting this integrated Eqn. (A.4) back into the von-Neumann equation (A.3), 
one obtains an integral form of the von-Neumann equation which reads 

p{t) = -t[Hj{t),pm - fds[Hi{t), [Hj{s),p{s)]]. (A.5) 

Jo 

For an observable O, the expectation value at a given time can be expressed as, 
= Trp(t)0. Therefore, using the integral form of the von-Neumann equation, 
the rate of change of O can be written as, 

dt{d), = - f ds T,{d[Hj{t), [H:{slp{s)]]}. (A.6) 
Jo 
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Here we have assumed that the initial density matrix p(0) commutes with the 
observable O. Then, the first term in Eqn. (A.5) yields a zero value when performing 
the trace. Expanding the double commutators and using the cyclic property of the 
trace to interchange density matrix to the rightmost side of each term, we find 

9t{d),^ - f dsTr{dHi(t)Hi(s)p(s) + Hi(s)Hi{t)dp{s) 
Jo 

- Hj{s)dHi{t)p{s) - Hi{t)dHi{s)p{s)}. (A.7) 
For — Pn this equation coincides with (6). 
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